GL-TR-89-017 1 


A  Multichannel  Processing  Approach  to  Real  Time 
Network  Detection,  Phase  Association  and 
Threshold  Monitoring 


Frode  Ringdal 
Tormod  Kvaerna 


AD-A213  940 


NTNF/NORSAR 
Post  Box  51 

N-2007  Kjeller,  NORWAY 


31  May  1989 


Scientific  Report  No.  1 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


GEOPHYSICS  LABORATORY 
AIR  FORCE  SYSTEMS  COMMAND 
UNITED  STATES  AIR  FORCE 

HANSCOM  AIR  FORCE  BASE,  MASSACHUSETTS  01731-5000 


SPONSORED  BY 

Defense  Advanced  Research  Projects  Agency 
Nuclear  Monitoring  Research  Office 
ARPA  ORDER  NO.  5307 

MONITORED  BY 
Geophysics  Laboratory 
Contract  No.  F49620-89-C-0038 


The  views  and  conclusions  contained  in  this  document  are  those  of  the  authors 
and  should  not  be  interpreted  as  representing  the  official  policies,  either 
expressed  or  implied,  of  the  Defense  Advanced  Research  Projects  Agency  or  the 
U.S.  Government. 


This  technical  report  has  been  reviewed  and  is  approved  for  publication. 


?S  F.  LEWKOWICZ 
rract  Manager 

io l  id  Earth  Geophysics  Branch 
Earth  Sciences  Division 


JAMES]  F.  LEWKOWICZ 
Bxgrigh  Chief 

Solid  Earth  Geophysics  Branch 
Earth  Sciences  Division 


FOR  THE  COMMANDER 


DONALD  H.  ECKHARDT,  Director 
Earth  Sciences  Division 


This  report  has  been  reviewed  by  the  ESD  Public  Affairs  Office  (PA)  and  is 
releasable  to  the  National  Technical  Information  Service  (NTIS). 

Qualified  requestors  may  obtain  additional  copies  from  the  Defense  Technical 
Information  Center.  All  others  should  apply  to  the  National  Technical 
Information  Service. 

If  your  address  has  changed,  or  if  you  wish  to  be  removed  from  the  mailing 
list,  or  if  the  addressee  is  no  longer  employed  by  your  organization,  please 
notify  AFGL/DAA,  Hanscom  AFB,  MA  01751-5000.  This  will  assist  us  in  main¬ 
taining  a  current  mailing  list. 


Do  not  return  copies  of  this  report  unless  contractual  obligations  or  notices 
on  a  specific  document  requires  that  it  be  returned. 


Unclassif  ied 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  _ 


REPORT  DOCUMENTATION  PAGE 


lb  RESTRICTIVE  MARKINGS 


Form  Approved 
OMB  No.  0704-0188 


Id  REPORT  SECURITY  CLASSIFICATION 

Unclass  if ied 


2a  SECURITY  classification  authority 


2b  DECLASSIFICATION  /  DOWNGRADING  SCHEDULE 


4  PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 


3  DISTRIBUTION /AVAILABILITY  OF  REPORT 
Approved  for  public  release; 
Distribution  unlimited. 


5  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 
GL-TR-89-0 171 


6a  NAME  OF  PERFORMING  ORGANIZATION 


6c.  ADDRESS  (City.  State,  and  ZIP  Code) 
Post  Box  51 

N-2007  Kj oiler,  Norway 


8c.  ADDRESS  (City,  State,  and  ZIP  Code) 


6b  OFFICE  SYMBOL  7a  NAME  OF  MONITORING  ORGANIZATION 
(If  applicable) 

Gcophys ics  Labor a  Tory 


7b  ADDRESS  (City,  State,  and  HP  Code) 

Hanscoin  Air  Force  Base 
Massachusetts  01731-5000 


9  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 
Contract  No.  F49620 - 89 - C - 0038 


10.  SOURCE  OF  FUNDING  NUMBERS 


1400  Wilson  Blvd. 

Arlington.  VA  22209-2308  _ 


1 1  TITLE  (Include  Security  Classification) 

A  Multichannel  Processing  Approach  to  Real  Time  Network  Detection,  Phase 
Association  and  Threshold  Monitoring  _ _ 


12  PERSONAL  AUTHOR(S) 

Frode  Ringdal  and  Tormod  Kvaerna 


13a.  TYPE  OF  REPORT 


13b  TIME  COVERED 

14  DATE  OF  REPORT  (Year,  Month,  Day) 

FROM  gq/02./aiTO 

fxvi%vwMH\ 

1989  May  31 

COSATI  CODES 


18  SUBJECT  TERMS  (Continue  on  reverse  if  necessary  and  identify  by  block  number) 


SUB-GROUP 


•/Network  detection,  phase  association, 

threshold  monitoring,  regional  arrays  • 


19  ABSTRACT  ( Continue  on  reverse  if  necessary  and  identify  by  block  number) 

This  paper  presents  an  approach  to  multichannel  processing  of  data  recorded  by  a  net¬ 
work  of  stations  which  represents  an  extension  of  the  delay-and-sum  beamforming  tradi¬ 
tionally  applied  in  array  processing.  A  number  of  geographical  beam- steering  points  are 
defined,  and  for  each  beam  a  set  of  time-aligned  traces  derived  from  the  network  sta¬ 
tions  are  processed  together  so  as  to  extract,  for  each  step  in  time,  a  set  of  features 
corresponding  to  that  particular  beam.  Applying  this  approach  to  the  network  of  the 
three  regional  arrays  NORESS,  ARCESS  and  FINESA  in  Fennoscandia ,  we  demonstrate  its 
usefulness  for  associating  regional  phases  detected  at  individual  arrays  and  origi¬ 
nating  from  the  same  event.  We  also  give  an  example  of  application  addressing  the 
problem  of  continuously  monitoring  the  seismic  noise  field.  In  this  regard,  we  show 
that  one  can  obtain,  at  a  given  confidence  level,  a  continuous  assessment  of  the  upper 
limit  of  magnitudes  of  seismic  events  that  would  go  undetected  by  such  a  network. 


20  DISTRIBUTION /  AVAILABILITY  OF  ABSTRACT 
Cl  UNCLASSl F lED/UNL IMI  TE  D  □  SAME  AS  RPT 
ZZa  NAME  OF  RESPONSIBLE  iNDiVIOUAi 

James  Lewkowicz 


□  DTIC  USERS 


21  ABSTRACT  SECURITY  CLASSIFICATION 

_ IINr.lAS.SI.FJ.Efl - __ - 

22b  TELEPHONE  (Include  Area  Code)  22c  OFFICE  SYMBOL 

(617)  37  /  -  3028  CL/LWH 


DO  Fo/m  1473.  JUN  86 


Previous  editions  are  obsolete 


SECURITY  Cl  ASSIF ICATION  OF  THIS  PAGE 


Pro  face 
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within  a  wide  range  of  subjects  relevant  to  seismic  monitoring.  The  emphasis  of  the 
research  program  is  on  developing  and  assessing  methods  for  processing  of  data 
recorded  by  networks  of  small-aperture  arrays  and  .'{-component  stations,  for  events 
both  at  regional  and  tcleseismic  distances.  In  addition,  more  general  seismological 
research  topics  are  addressed. 

I'ach  quarterly  technical  report  under  this  contract  will  present  one  or  several 
separate  investigations  addressing  specific  problems  within  hie  scope  of  the  state¬ 
ment  of  work.  Summaries  of  the  research  t  hurts  within  the  program  as  a  whole  will 
be  given  in  annual  reports. 

This  Scientific  Ilcpcrt  No.  1  pitociiw  a  manuscript  entitled  “A  multichannel 
processing  approach  to  real  time  network  detection,  phase  association  and  threshold 
monitoring",  by  Frode  Ringdal  and  Tormod  Kvanrna. 
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A  multichannel  processing  approach 
to  real  time  network  detection,  phase 
association  and  threshold  monitoring 


FRODE  RINGDAL  and  TOKMUD  KV/ERNA 
NTNF/NORSAR,  Post  Box  51,  N-2007  Kjoller,  Norway 


Short  Title:  Multichannel  processing  of  network  data 


Abstract 


This  paper  presents  an  approach  to  multichannel  processing  of  data  recorded  by  a 
network  of  stations  which  represents  an  extension  of  the  delay-and-sum  beamforming 
traditionally  applied  in  array  processing.  A  number  of  geographical  beam-steering 
points  are  defined,  and  for  each  beam  a  set  of  time-aligned  traces  derived  from 
the  network  stations  are  processed  together  so  to  extract,  for  each  step  in  time, 
a  set  of  features  corresponding  to  that  particular  beam.  Applying  this  approach 
to  the  network  of  the  three  regional  arrays  NOR  ESS,  ARCESS  and  FINESA  in 
Fennoscandia,  we  demonstrate  its  usefulness  for  associating  regional  phases  detected 
at  individual  arrays  and  originating  from  the  same  event.  We  also  give  an  example 
of  application  addressing  the  problem  of  continuously  monitoring  the  seismic  noise 
field.  In  this  regard,  we  show  that  one  can  obtain,  at  a  given  confidence  level,  a 
continuous  assessment  of  the  upper  limit  of  magnitudes  of  seismic  events  that  would 
go  undetected  by  such  a  network. 

Introduction 


In  the  processing  of  seismic  network  data,  individual  phase  detections  corre¬ 
sponding  to  the  same  seismic  event  must  be  properly  associated  and  grouped  to¬ 
gether.  This  is  today  usually  done  starting  with  an  initial  trial  epicenter  and  then 
applying  various  search  strategies  supplemented  by  combinational  techniques.  For 
teleseismic  monitoring  using  global  network  data,  such  techniques  are  well  estab¬ 
lished,  and  a  large  degree  of  automation  has  been  achieved  (Goncz,  1980;  Slunga, 
1980). 

In  recent  years,  the  subject  of  regional  monitoring  using  networks  comprising 
small-aperture  arrays  as  well  as  single  stations  has  attracted  increased  attention. 
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The  motivation  has  been  the  need  for  improving  the  capability  to  monitor  under¬ 
ground  nuclear  explosion  testing  down  to  very  low  magnitudes.  The  inclusion  of 
regional  phases  in  the  phase  association  procedure  leads  to  a  considerable  increase 
in  the  complexity  of  this  task,  and  much  effort  has  been  devoted  to  developing 
automated  procedures  for  handling  this  problem. 

Techniques  for  automatic  association  of  regional  seismic  phases  recorded  by  a 
single  array  were  iirst  established  in  connection  with  the  early  N011ESS  array  de¬ 
velopments  (Mykkeltveit  and  Itingdai,  1981;  Mykkeltveit  and  Bungum,  1981).  In 
later  developments,  knowledge-based  system  concepts  were  introduced,  and  the  al¬ 
gorithms  were  expanded  to  include  processing  of  data  from  a  network  of  regional 
arrays  and  single  stations  (Bache,  1987;  Baumgardt,  1987).  Bratt  and  Bache  (19SS) 
developed  an  automated  procedure  for  locating  regional  seismic  events  recorded  by 
such  a  network,  incorporating  arrival  times  and  azimuth  information. 

The  emphasis  of  the  phase  association  methods  developed  so  far  has  been  to 
build  on  the  techniques  previously  applied  successfully  for  global  teleseismic  net¬ 
works,  i.e.,  to  associate  individual  phase  detections  at  network  stations  using  combi¬ 
national  techniques.  Multichannel  processing  methods  have  received  little  attention, 
mainly  because  the  signal  coherency  across  extended  networks  is  too  low  to  make 
conventional  beamformiug  useful.  Incoherent  beamforming  (Ringdal  et  al,  1972; 
Husebye  et  al,  1972)  has  been  suggested  as  a  possible  alternative  approach.  How¬ 
ever,  whereas  this  method  ic  effective  for  networks  of  limited  aperture,  its  practical 
usefulness  for  larger  networks  remains  uncertain. 

This  paper  presents  a  multichannel  processing  approach  for  network  data  that 
we  have  termed  ‘‘generalized  beamforming”.  In  a  sense,  it  represents  not  one  single 
technique,  but  rather  a  framework  for  processing  such  data.  Applying  a  procedure 
similar  to  the  conventional  delay-and-sum  beamforming  used  in  array  processing,  we 
define  a  number  of  geographical  beam-steering  points  covering  the  region  of  interest. 
For  each  beam  point,  we  obtain  a  set  of  time-aligned  input  traces,  from  which  a  set  of 
beam  features  are  continuously  extracted.  We  show  how  this  approach,  in  case  of  a 
regional  network,  can  be  used  to  associate  individual  phase  detections  corresponding 
to  the  same  event.  We  also  give  an  example  of  application  addressing  the  problem 
of  continuously  monitoring  the  seismic  noise  level,  for  the  purpose  of  obtaining  a 
quantitative  assessment  of  the  upper  limit  of  magnitudes  of  seismic  events  that 
would  go  undetected  by  such  a  network. 

Method 

Iri  this  section,  we  outline  the  overall  approach  to  the  problems  addressed  in 
this  paper.  We  do  not  go  into  detail  at  this  stage,  but  refer  instead  to  subsequent 
sections  discussing  how  the  general  framework  can  be  used  in  practical  application 
to  various  seismic  monitoring  problems. 

Let  us  assume  that  a  network  of  N  seismic  stations  is  available  for  monitoring 
a  specified  geographical  region.  For  simplicity  of  presentation,  we  will  assume  that 


these  are  all  array  stations,  able  to  provide  phase  velocity  and  azimuth  information 
for  detected  signals.  Extension  to  the  single-station  case  is  straightforward,  although 
the  quality  of  the  results  will  naturally  be  reduced  compared  to  when  array  stations 
are  available. 

We  first  subdivide  the  region  to  bo  monitored  by  introducing  a  grid  of  J  ge¬ 
ographical  aiming  points.  Each  such  point  then  corresponds  to  a  beam  location, 
so  that  the  beam  set  covers  the  entire  region  with  a  predefined  spacing  of  the  grid 
points. 

For  each  station  in  the  network,  we  assume  tha*  the  recorded  data  are  processed 
separately,  using  conventional  algorithms.  Specifically,  let  us  denote  by  .s, , ^.( 7’) 
the  beam  at  the  i'th  array  (i  =  1,2,...  A’),  steered  toward  the  j’tli  beam  point 
(j  =  1,2,. ..J)  and  corresponding  to  the  k'  th  seismic  phase  [k  =  1,2, . . . ,  A\j). 
Here,  KtJ  denotes  the  number  of  phases  that  might  be  of  interest  for  the  particular 
station-beam  combination.  We  assume  that  these  traces  are  subjected  to 

standard  automatic  detection  processing,  thus  producing  lists  of  signal  onset  times, 
phase  velocity/azimuth  estimates  and  other  detection  parameters,  as  well  as  noise 
level  estimates  during  periods  of  non-detection. 

In  analogy  with  conventional  array  beamforming,  the  procedure  is  now  to  “steer” 
the  network  toward  each  beam  location  and  process  each  beam  individually  by 
continuously  extracting  relevant  features  for  tha!  beam  as  a  function  of  time. 

The  network  beamsteering  is  done  by  computing  a  se*  of  time  delays  r,y/t,  with 
indices  defined  as  before,  for  all  combinations  of  beams,  stations  and  phases.  Stan¬ 
dard  travel-time  tables  are  used  in  these  computations.  Thus,  for  the  j’th  beam,  we 
obtain  a  set  of  time-aligned  channels: 

-Sj(T)  =  {s,jk(T  +  rt}k)}  k  =  1....  luy,  (1) 

Here,  $j{T)  can  be  viewed  as  a  signal  vector  representing  the  individual  sta¬ 
tion/phase  observations  corresponding  to  a  hypothetical  event  with  origin  time  T 
and  located  at  the  aiming  point  of  the  j’th  beam. 

Given  this  time  alignment  of  the  input  traces,  we  now  can  proceed  in  several 
different  ways  to  extract,  as  a  function  of  time,  characteris-  tic  features  of  the  seismic 
field  as  seen  when  focusing  on  a  specific  beam  area.  One  such  approach  is  to  consider 
the  detector  outputs  on  each  individual  trace  and  combine  this  information.  In 
this  case,  the  network  beamforming  process,  for  a  given  beam,  at  time  T,  can  be 
described  as  looking  for  a  pattern  of  detections/non-detections  that  matches  the 
predicted  pattern  for  a  hypothetical  event  with  origin  time  T  and  location  within 
the  beam  region.  The  actual  beam  value  is  derived  from  probabilistic  considerations, 
and  in  essence  describes  how  well  the  observed  pattern  matches  the  prediction.  By 
moving  along  the  time  axis,  we  thus  obtain  a  beam  trace  that  can  be  subjected 
to  standard  threshold  algorithms  for  detection.  The  process  can  be  supplemented 
by  various  individual  "quality  of  fit”  measures  calculated  at  each  time  point.  An 
example  of  this  type  of  approach  for  the  purpose  of  regional  phase  association  is 
given  in  the  following. 
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The  generalized  beamforming  approach  also  provides  a  convenient  tool  to  con¬ 
tinuously  assess  the  seismic  noise  field  associated  with  a  given  beam.  An  application 
of  particular  interest  in  a  monitoring  situation  would  be  to  calculate,  at  each  step  in 
time,  upper  confidence  limits  for  the  magnitude  of  possible  non-detected  events  for 
each  beam.  This  would  be  useful  to  obtain  a  realistic  assessment  of  actual  network 
detection  capabilities,  at  any  given  point  in  time.  The  paper  presents  an  example 
of  practical  application  of  this  approach. 


Regional  phase  association 


The  method  has  been  applied  to  a  data  base  comprising  24  hours  of  recordings 
from  the  regional  arrays  NORESS,  ARCESS  and  FINESA  (Mykkeltveit  et  al,  1987; 
Korhonen  et.  ai.  1987;  Kvaerna,  1989),  with  a  beam  deployment  covering  Fennoscan- 
dia  and  adjacent  areas. 

A  RONAPP-typc  detector  (Mykkeltveit  and  Bungum,  1984)  was  first  applied  to 
each  array  individually,  using  the  broad-band  F-K  method  (Kvaerna  and  Doornbos, 
1986)  to  obtain  phase  velocity  and  azimuth  for  each  detected  phase.  The  resulting 
detection  lists  then  provided  the  input  to  the  network  processor. 

The  beam  grid  used  for  network  processing  is  shown  in  Figure  1,  and  comprises 
altogether  121  aiming  points,  approximately  equally  spaced.  Typical  distance  be¬ 
tween  aiming  points  is  150  km. 

In  the  network  beamforming  process,  a  simple  model  of  assigning  0/1  proba¬ 
bilities  to  individual  phases  at  each  station  was  used.  We  required  that  estimated 
phase  velocities,  azimuth,  dominant  frequency  and  arrival  times  fall  within  prede¬ 
fined  ranges  for  a  phase  detection  to  be  accepted  for  a  given  beam.  These  tolerance 
ranges  are  specified  in  Table  1.  Note  in  particular  that  only  very  general  criteria  are 
applied,  and  we  have  made  no  attempt  to  optimize  performance  by  regionalization. 

With  this  simplified  model,  the  network  beamforming  process  in  practice  was 
reduced  to,  for  each  beam  and  each  time  T,  counting  the  number  of  phase  matches 
for  a  hypothetical  event  located  in  the  beam  region  and  having  origin  time  T.  The 
detection  threshold  was  set  equal  to  2.  Thus,  all  occurrences  of  two  or  m  re  matching 
phase  detections  were  flagged  as  potential  events.  A  typical  beam  ..race  is  shown  in 
Figure  2. 

In  analogy  with  conventional  array  processing,  the  beamforming  procedure  oc¬ 
casionally  produces  side  lobe  detections,  thus  resulting  in  several  different  beam 
detections  for  a  given  event.  A  grouping/  reduction  process  is  therefore  required. 

The  grouping  procedure  applied  in  our  case  consisted  of  successively  linking 
together  entries  in  the  overall  beam  detection  list.  This  was  dune  in  such  a  way 
that  a  new  entry  would  be  linked  if  it  had  at  least  one  individual  phase  detection  in 
common  with  a  previous  entry  in  the  group.  The  maximum  allowable  duration  of 
a  group  was  set  to  10  minutes  (in  practice,  the  longest  duration  was  7  minutes  for 
this  data  set).  In  order  to  resolve  obvious  multiple  events,  groups  were  solit  up  if 
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two  P-detections  from  the  same  array  occurred  with  more  than  30  seconds  arrival 
time  difference. 

The  results  are  summarized  in  Table  2.  It  is  important  to  note  that  the  total 
of  91  groups  comprise  all  possible  events  that  could  be  associated,  given  the  station 
detection  lists.  Also,  a  scrutiny  of  the  data  shows  that  only  3  of  these  groups 
contain  multiple  events,  all  of  these  being  small  presumed  mining  explosions  seen 
by  one  array  only. 

Some  of  the  entries  in  Table  2,  c.g.,  those  generated  from  t  wo  secondary  phases, 
are  probably  questionable  seismic  events,  and  even  if  real,  may  be  impossible  to 
locate  accurately  without  access  to  additional  data.  An  upper  magnitude  limit 
could  be  estimated  for  such  events,  in  order  to  determine  whether  further  detailed 
analysis  is  desirable.  However,  the  large  majority  of  the  entries  appear  to  correspond 
to  real  seismic  events,  and  the  grouping  procedure  facilitates  the  subsequent  detailed 
analysis  of  the  associated  phases. 

The  network  beamforming  procedure  gives  an  initial  estimate  of  event  location 
by  selecting  the  ’’best  beam”  in  each  group.  'Phis  is  defined  as  the  beam  with  the 
greatest  number  of  associated  phase  detections,  and  if  equality,  the.  smallest  average 
time  residual  of  the  detected  phases.  Since  the  initial  beam  grid  is  very  coarse,  we 
applied  a  beampacking  algorithm  for  each  detection  group,  using  a  grid  spacing  of 
20  km  in  order  to  improve  the  location  estimate.  T  he  resulting  location  estimates 
for  the  data  set  are  displayed  in  Figure  3. 

Table  3  lists  the  results  of  the  automatic  procedure,  after  beam-  packing,  for 
those  events  for  which  independent  location  estimate's  based  on  local  network  data 
were  available.  We  note  that  the  respective  estimates  are  very  consistent  (median 
difference  40  km),  and  thus  the  beam  results  can  be  used  as  a  reasonable  first  esti¬ 
mate  of  event  location.  For  more  accurate  results,  available  techniques  for  accurate 
hypocenter  location,  e.g.,  the  TTAZLOC  procedure  (Bratt  and  Bache,  1988))  should 
be  used. 

Continuous  monitoring  of  upper  event  magnitude  limits 

As  a  second  application  of  the  generalized  beamforming  procedure,  we  now  ad¬ 
dress  the  problem  of  monitoring  the  noise  levels  on  each  bnam.  an^  "c“  this  mfor- 
mation  to  assess  the  size  of  events  that  might  go  undetected. 

In  formulating  the  approach,  we  consider  a  given  geographical  location,  and  a 
given  “origin  time”  of  a  hypothetical  event.  Assume  that  N  seismic  phases  are 
considered  (there  might  be  several  stations  and  several  phases  per  station). 

For  each  phase,  we  assume  that  we  have  an  estimate  5,  of  the  signal  (or  noise) 
level  at  the  predicted  arrival  time.  For  P-phases,  S ,  might  be  the  maximum  short 
term  average  (STA)  value  (1  second  integration  window)  within  +  5  seconds  of  the 
predicted  time.  For  Lg,  a  longer  STA  integration  window  (e.g.,  10  seconds)  might 
be  used,  and  its  maximum  might  be  selected  allowing  a  somewhat  greater  deviation 
from  the  predicted  arrival  time. 


YV  e  assume  that  the  network  has  been  calibrated  (or  alternatively  that  stan¬ 
dard  attenuation  values  are  available),  so  that  magnitude  correction  factors  (6.)  are 
available  for  all  phases.  Thus,  if  a  detectable  signal  is  present: 

m,  -  log(Si)  +  b,  (t=  1,2,  —  JV)  (2) 

Here.  mt  are  estimates  of  the  event  magnitude  m.  Statistically,  we  can  consider 
each  m,  as  sampled  from  a  normal  distribution  [in, a).  Based  on  NORSAR  experi¬ 
ence,  we  consider  a  standard  value  of  a  =  0.2  to  be  reasonable  for  a  small  epicentral 
area,  and  this  value  will  be  used  in  the  following. 

Let  us  now  assume  a  “noise  situation”,  i.e.,  that  there  are  no  phase  detections 
corresponding  to  events  at  the  given  location  for  the  given  origin  time. 

We  then  have  a  set  of  “noise”  observations  a,,  where 

",  =  log(St)  +  bt  (i  =  1,2 . V)  (3) 

If  a  hypothetical  event  of  magnitude  m  were  present,  it  would  have  phase  mag¬ 
nitudes  m,  normally  distributed  around  m.  We  know  that  for  e=>rh  phase, 

mt  <  a,  (i  =  1,2, ...  N)  (4) 

Following  a  proceduie  similar  to  that  of  Ringdal  (1976),  we  now  consider  the 
function: 

f{rti)  =  Prob(all  m,  <  a,  /  event  magnitude  m)  (5) 

For  each  phase,  we  obtain  probability  functions  /,(m)  and  Sf,(m)  as  follows: 

/,(m)  =  Frobfm,  <  a,/m)  =  1  -  $(— — — )  (t  =  1,2, . . .  A')  (6) 

r 

g,(m)  =  ?Tob(m,  >  a,/m)  =  (i  =  1 , 2, . . .  .V)  (7) 

where  $  is  the  standard  (0,1)  normal  distribution. 

Thus,  assuming  independence, 

N 

/(m)  =  I]  /i(m)  (8) 

•=i 


The  probability  g(m)  that  at  least  one  of  the  observed  noise  values  would  be 
exceeded  by  the  signals  of  a  hypothetical  event  of  magnitude  m,  then  becomes 

g(m)  =  1  -  f(m)  (9) 

As  illustrated  in  Figure  4,  the  90  per  cent  upper  limit  is  then  defined  as  the 
solution  of  the  equation 

gym)  =  0.90  (10) 


h 


It  is  impor' .  ii  U'  interpret  the  90  per  rent  limit  d*  -lined  above  in  the  proper  way. 
Thus,  it  s>.  el  not  be  considered  as  a  9b  per  cent  network  detection  threshold  since 
we  huv  ...ade  no  allowance  for  a  signal-to-noise  ratio  which  would  be  required  in 
order  to  delect  an  event,  given  the  noise  levels.  Rather,  the  computed  level  is  tied  to 
the  actually  observed  noise  values,  and  to  the  fact  that  any  hypothetical  signal  must 
In'  below  these  values.  Our  90  per  cent  limit  represents  the  largest  magnitude  of  a 
possible  hidden  event,  in  the  sense  that  above  ties  limit,  there  is  at  least  a  90  per 
cent  probability  that  one  or  more  of  the  observed  noise  values  would  be  exceeded 
b\  the  signals  of  such  an  event. 

As  an  application  (,f  the  method,  we  selected  an  area  as  shown  in  Figure  5 
situated  at  similar  distance  from  the  tluoe  airay-.  for  each  of  the  three  arrays, 
one  l’n  beam  and  one  Lg  beam  were  .steer'd  to  this  location.  The  beam  traces 
were  filtered  using  the  frequency  bands  3  -5  11/  (I'n>  and  2  -4  Hz  (Lg).  Magnitude 
calibration  values  (/>,)  were  obtained  i>\  processing  previously  recorded  events  of 
known  magnitude  (M/  )  and  at  similar  distance  ranges,  and  then  determining  b, 
vaiti'-s  independently  for  Pn  and  l,g. 

Based  on  these  input  traces  from  the  three  arrays,  a  network  beam  was  then 
formed,  using  time  delays  for  eaeh  phase  that  corresponded  to  the  given  location. 
Arrival  time  tolerances  wore  set  to  ±  5  seconds  for  Pn  and  ±  10  seconds  for  Lg. 
This  is  roughly  consistent  with  a  beam  radius  of  50  km  as  shown  on  the  figure.  STA 
integration  windows  were  set  to  1  second  for  Pn  and  10  seconds  for  Lg.  The  values  of 
Si  in  eq.  (2)  were  obtained  as  the  maximum  STA  values  within  the  respective  arrival 
time  tolerances,  using  the  mid-point  of  the  integration  interval  as  time  reference. 

We  chose  to  analyze  a  3  1/2  hour  interval  during  which  four  regional  seismic 
events  of  M/,  >  2.0  were  reported  in  the  Helsinki  bulletin.  These  events  were  all 
located  outside  the  beam  region  to  be  studied,  and  one  of  our  aims  was  to  investigate 
how  interfering  signals  from  these  events  would  infiimnco  the  monitoring  capability 
for  the  chosen  beam  region. 

Figure  6  shows,  for  the  beam  region  considered,  the  computed  90  per  cent  upper 
magnitude  limits,  plotted  as  a  function  of  time.  In  this  figure,  only  the  Pn  phase  has 
been  used,  and  the  three  arrays  are  shown  individually  and  in  combination  (bottom 
trace). 

It  is  clear  from  Figure  6  that  when  considering  individual  arrays  only,  there  are 
several  possible  time  intervals  when  relatively  large  events  (Ml  ~  2. 0-3.0)  located  in 
the  beam  area  might  go  undetected  because  of  signals  from  interfering  events.  How¬ 
ever,  when  the  Pn  phases  are  combined,  these  instances  necur  much  more  seldom. 

Figure  7  shows  a  similar  plot,  but  this  time  including  both  the  Pn  and  the  Lg 
phase  for  each  array.  Even  on  an  individual  array  basis,  this  causes  substantial 
reduction  in  the  upper  magnitude  limits.  For  the  combined  plot  (bottom  trace  of 
Figure  7),  which  takes  into  account  all  6  Pn  and  Lg  phases  from  the  three  arrays,  we 
see  that  the  upper  limit  is  well  below  M/,  =  2.0  for  the  entire  time  interval.  Thus, 
we  may  conclude  that,  at  the  specified  level  of  lonfmence,  no  event  of  M/,  =  2.0  or 
higher  occurred  in  the  beam  region  during  the  time  period  considered. 


Discussion 


With  regard  to  phase  association,  the  generalized  beamforming  approach  pre¬ 
sented  in  tins  paper  provides  an  effective  method  to  group  all  combinations  of  in¬ 
dividual  phase  detections  that  could  possibly  correspond  to  the  same  seismic  event. 
At  the  same  time,  preliminary  estimates  of  epicenter  and  origin  time  are  obtained. 

The  primary  importance  of  this  would  be  to  obtain  a  starting  point  for  subse¬ 
quent  detailed  interactive  analysis  aimed  at  precise  determination  of  source  param¬ 
eters.  In  particular,  expert  system  approaches  (either  script-based  or  rule-based) 
could  be  invoked  at  this  stage.  The  advantage  of  applying  tire  generalized  beam- 
forming  as  the  first  step  is  to  reduce  the  amount  of  combinational  processing  that 
would  be  necessary  otherwise.  It  is  here  noteworthy  that  the  processing  load  when 
applying  beamforming  increases  in  a  linear  fashion  when  the  number  of  individ¬ 
ual  phase  detections  increase,  whereas  combinational  possibilities  tend  to  increase 
exponentially.  While  wo  have  in  this  paper  used  only  a  three-array  network,  the 
extension  to  larger  networks  is  clearly  straightforward. 

The  application  of  the  method  to  provide  continuous  monitoring  of  upper  magni¬ 
tude  limits  at  specified  beam  locations  provides  a  useful  supplement  to  standard  sta¬ 
tistical  network  capability  studies  (e.g.,  Wirth,  1977;  Ringdal,  1986).  In  particular, 
this  application  would  give  a  way  to  assess  the  possible  magnitude  of  non-detected 
events  during  the  coda  of  large  earthquakes.  In  such  situations,  it  would  be  appro¬ 
priate  to  use  global  network  data  and  include  as  many  relevant  phases  as  possible 
for  each  network  station.  For  example,  while  an  expected  P  phase  at  a  given  station 
may  be  obscured  by  the  earthquake  coda,  later  phases  such  as  PcP  or  PP  may  be 
less  influenced,  and  the  noise  level  at  their  respective  expected  arrival  times  would 
therefore  provide  important  information  as  to  the  size  of  possible  undetected  events. 

As  a  final  comment,  we  note  that  the  approach  presented  here  to  upper  limit 
magnitude  calculation  could  be  applied  to  extend  the  utility  of  various  discriminants, 
such  as  M,:mt.  For  small  explosions,  surface  waves  frequently  are  too  weak  to  be 
observed  at  any  station  of  the  recording  network.  Obtaining  reliable  upper  bound 
on  Mi  in  such  cases  would  expand  the  range  of  usefulness  of  this  discriminant.  In 
practice,  an  ’’upper  bound”  for  single-station  measurements  has  often  bre  ,  given  as 
the  ’’noise  magnitude”  at  that  station,  i.e.,  the  M,  value  that  cor r<- .ponds  to  the 
actually  observed  noise  level  at  the  expected  time  of  Rayleigh  wave  arrival.  The 
proposed  procedure  will  include  this  as  a  special  case  of  a  more  general  network 
fori.  Ration. 

In  future  studies,  we  plan  to  investigate  the  application  of  more  sophisticated 
probabilistic  models  in  the  generation  of  beam  traces  and  the  continuous  extraction 
of  features  associated  with  the  individual  beams.  Application  to  larger  networks, 
including  teleseismic  monitoring  using  global  network  data,  will  also  be  considered. 
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Table  Captions 


Table  1.  Acceptance  limits  for  parameters  used  in  the  network  beamforming  pro¬ 
cess,  as  applied  in  the  example  described  in  the  text. 

Table  2.  Phase  groups  associated  by  the  network  beamforming  procedure  for  a 
24-hour  interval,  based  on  NORESS/ARCESS/FINESA  detection  lists. 

Table  3.  Location  estimates  obtained  automatically  from  the  beampacking  proce¬ 
dure  compared  with  independent  network  locations  from  the  Helsinki  and  Bergen 
bulletins.  Note  the  good  consistency,  especially  for  events  with  more  than  one  de¬ 
tecting  array. 
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Figure  Captions 


Fig.  1.  Beam  grid  used  in  the  generalized  beamforrning  procedure  for  the  purpose 
of  associating  regional  phases  from  NORESS,  ARCESS  and  E1NESA.  The  location 
of  the  three  arrays  is  shown  on  the  map. 

Fig.  2.  Illustration  of  the  beamforming  procedure,  using  0/1  weights  for  individual 
phase  detections  as  described  in  the  text.  For  each  of  the  three  arrays,  detection 
time  traces  for  Pn  and  Lg  are  displayed  in  the  form  of  step  functions.  A  value  of  1 
for  these  functions  indicates  that  an  acceptable  detection  (with  regard  to  azimuth, 
phase  velocity,  etc.)  has  occurred  within  a  specified  time  window.  To  illustrate  the 
beam  delays,  an  arbitrary  origin  timeT  has  been  marked  on  the  beam  trace,  and  the 
predicted  arrival  times  corresponding  to  a  hypothetical  event  at  the  beam  location 
with  origin  time  T  are  marked  as  arrows.  The  network  beam  (top  trace)  is  a  sum  of 
the  time-aligned  individual  traces,  including  the  Pg  and  Sn  phases,  which  are  not 
shown  on  the  plot. 

Fig.  3.  Event  location  results,  after  beampacking,  for  the  phase  groups  associated 
by  the  network  beamforming  algorithm.  The  location  of  the  three  arrays  NORESS, 
ARCESS,  FINESA  is  also  shown. 

Fig.  4.  Illustration  of  the  procedure  for  calculating  upper  magnitude  limits.  Each 
network  station  gives  rise  to  a  probability  distribution  gt(M)  as  described  in  the 
text.  'The  dottetTcurve,  g(M),  represents  the  probability,  given  event  magnitude  A/, 
that  the  signal  from  a  hypothetical  event  would  exceed  the  actually  observed  noise 
level  at  at  least  one  station. 

Fig.  5.  Location  of  the  beam  area  used  in  the  example  of  continuous  monitor¬ 
ing  of  upper  magnitude  limits  on  non-detected  events.  The  area  covers  a  circle  of 
approximately  50  km  radius,  and  is  situated  at  similar  distances  from  the  three 
arrays. 

Fig.  6.  Results  from  the  continuous  threshold  monitoring  of  the  area  shown  in 
Figure  3  for  a  3  1/2  hour  period,  using  Pn  phases  only.  The  top  three  traces  show, 
for  each  array,  the  largest  magnitude  of  a  possible  non-detected  event  (confidence 
90  per  cent)  as  a  function  of  time.  The  bottom  trace  shows  the  result  of  combining 
the  observations  from  all  three  arrays  (Pn  phase  only)  as  described  in  the  text. 

Fig.  7.  Same  as  Figure  6,  but  using  both  the  Pn  and  Lg  phases  for  the  upper 
magnitude  limit  calculations.  Comparing  with  Figure  6,  we  note  that  this  serves  to 
lower  the  thresholds,  both  for  each  individual  array  (top  three  traces)  and  for  the 
combined  results  (bottom  trace). 
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Table  1 


Phase  Type 

Pn  Pg  Sn  Lg  Rg 


Distance  interval^-) 
(kra)  for  which  a 
phase  is  accepted 

160-3000 

0-600 

160-3000 

0-2000 

0-A00 

Maximum  allowable 
deviation  from 
predicted  arrival 
time  (s) 

15 

20 

30 

35 

40 

Maximum  allowable 
azimuth  deviation 
(degrees) 

20 

.  20 

20 

20 

20 

Acceptance  limits 
for  apparent  phase 
velocity  (km/s) 

5.8-14 

5.8-10 

3. 2-5. 8 

3. 0-5.0 

2. 5-3. 7 

Acceptance  limits^ ) 
for  dominant 

0.5-20 

0.5-20 

0.5-20 

0.5-20 

0.5-20 

frequency  (Hz) 


"  NORESS,  the  Rg  phase  is'  not  included  in  the  phase  table 

For  FINESA,  a  lower  frequency  limit  of  0.9  Hz  is  used  for  all  phases. 


1 A 


Table  2 


Number  of  phases  for  best  beam 
in  each  group 


Number  of  phase 

groups : 

2 

3 

4 

5 

6 

7 

8 

NORESS  only 

18 

13 

4 

1 

0 

0 

0 

0 

ARCESS  only 

34 

19 

10 

4 

1 

0 

0 

0 

FINESA  only 

14 

13 

1 

0 

0 

0 

0 

0 

Two  arrays 

17 

9 

4 

3 

0 

1 

0 

0 

Three  arrays 

8 

0 

0 

2 

0 

1 

3 

2 

Totals 

91 

54 

19 

10 

1 

2 

3 

2 
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Table  3 
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LATITUDE  (DEG  N) 


Beam  location  57.0N  10.5E 


077/08:20:00.000 


PROBABILITY 
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Figure  4 


Probability  that  signal  exceeds  observed  noise  ,  given  M 
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Figure  5 
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Washington,  D.C.  20451 
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